Efficient Simulations of Interstellar Gas-Grain Chemistry 

Using Moment Equations 
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ABSTRACT 



O . 

O ' Networks of reactions on dust grain surfaces play a crucial role in the chem- 

istry of interstellar clouds, leading to the formation of molecular hydrogen in 
diffuse clouds as well as various organic molecules in dense molecular clouds. 
(-( ■ Due to the sub-micron size of the grains and the low flux, the population of re- 

0-^• active species per grain may be very small and strongly fluctuating. Under these 

O ' conditions rate equations fail and the simulation of surface-reaction networks 

■ requires stochastic methods such as the master equation. However, the master 

^ ■ equation becomes infeasible for complex networks because the number of equa- 

tions proliferates exponentially. Here we introduce a method based on moment 
>- ■ equations for the simulation of reaction networks on small grains. The number 

of equations is reduced to just one equation per reactive specie and one equation 



CN ' per reaction. Nevertheless, the method provides accurate results, which are in 

^ ■ excellent agreement with the master equation. The method is demonstrated for 



the methanol network which has been recently shown to be of crucial importance. 

Subject headings: dust — ISM; abundances — ISM; molecules — ISM; clouds — 
ISM; molecular processes 



1. Introduction 

Chemical networks in iri t erstellar cloud s consist of gas-phase and grain-surface reactions 



( iHartquist fc Williams 19951 : iTielens 20051 ). Reactions that take place on dust grains in 



elude the formation of molecular hydrogen (ICould fc Salpeter 19631 : iHoUenbach et al. 19711 ) 



as well as reaction networks producing ice mantles and various organic molecules. Unlike 
gas phase reactions in cold clouds that mainly produce unsaturated molecules, surface pro- 
cesses are dominated by hydrogen-addition reactions that result in saturated, hydrogen-rich 
molecules, such as H2CO, CH3OH, NH3 and CH4. In particular, rec ent experiments show 



that methanol cannot be efficiently produced by gas phase reactions (ICeppert et al. 20061 ). 
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On the other hand, the re are indications that it can be efficiently produced on ice-coated 
grains (IWatanabe 20051 ). Therefore, the abihty to perform simulations of the production o f 
methanol and other complex molecules on grains is of great import ance (IGarrod et al. 20061) . 
Unlike gas-phase reac tions, simulated using rate equation models (jPickles fc Williams 19771 : 

Hasegawa et al. 19921). grain-surface reactions r equire stochastic methods such as t he master 

equation (IBiham et al. 200ll : lGreen et al. 20011 ). or Monte Carlo (MC) simulations (jCharnley 20011 ) 
This is due to the fact that under interstellar conditions, of extremely low gas density 
and sub-micron grain sizes, surface reaction rates are dominated by fluctuations which 
cannot be accounted for by rate equati ons (ITielens fc Hagen 19821 : ICharnley et al. 19971 : 
Caselli et al. 19981 : IShalabiea et al. 19981 ). A significant advantage of the master equation 
over MC simulations is that it consists of differential equations, which can be easily cou- 
pled to the rate equations of gas-phase chemistry. Furthermore, unlike MC simulations 
that require the accumulation of statistical information over long times, the master equa- 
tion provides the probability distribution from which the reaction rates can be obtained 
directly. However, the number of equations increases exponentially with the number of reac 



five species, making the si mulation of complex networks infeasible ( IStantcheva et al. 2002 



Stantcheva fc Herbst 20031 ). The recently proposed multi-plane method dramatically reduces 
the number of equations, by breaking the network into a set of fully connected sub-networks 
( iLipshtat fc Biham 20041 ). enabling the simulation of more complex networks. However, the 
construction of the multi-plane equations for large networks turns out to be difficult. 

In this Letter we introduce a method based on moment equations which exhibits crucial 
advantages over the multi-plane method. The number of equations is further reduced to the 
smallest possible set of stochastic equations, including one equation for the population size of 
each reactive specie (represented by a first moment) and one equation for each reaction rate 
(represented by a second moment). Thus, for typical sparse networks the complexity of the 
stochastic simulation becomes comparable to that of the rate equations. Unlike the master 
equation (and the multi-plane method) there is no need to adjust the cutoffs - the same set 
of equations applies under all physical conditions. Unlike the multi-plane equations, the mo- 
ment equations are linear and for steady state conditions can be easily solved using algebraic 
methods. Moreover, for any given network the moment eq uations can be easily c onstructed 
using a diagrammatic approach, which can be automated (IBarzel fc Biham 20071 ). 



2. The Method 



To demonstrate the method we consider a simple network, sho wn in Fig. [H that 
involves three reactive species: H and O atoms and OH molecules (ICaselli et al. 19981 : 
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Shalabiea et al. 19981 : IStantcheva et al. 20021 ) . For simplicity we denote the reactive species 



by Xi = H, X2 = O, X3 = OH, and the non-reactive product species by X4 = H2, 
-^5 = O2, = H2O. The reactions that take place in this network include H + O ^ 
OH (Xi + X2 ^ X3), H + H ^ H2 (Xi + Xi ^ X4), O + O O2 (X2 + X2 X5), and 
H + OH ^ H2O (Xi + X3 -> Xe). 

Consider a spherical grain of diameter d, exposed to fluxes of H and O atoms and OH 
molecules. The cross-section of the grain is a = tkP/A and its surface area is ird?. The 
density of adsorption sites on the surface is denoted by s (sites cm~^). Thus, the number 
of adsorption sites on the grain is S* = ttc/^s. The desorption rates of atomic and molecular 
species from the grain are given hj Wi = u ■ exp[—Ei{i)/kBT], where u is the attempt rate 
(standardly taken to be 10^^ s"^), Ei{i) is the activation energy for desorption of specie Xj 
and T (K) is the grain temperature. The hopping rate of adsorbed atoms between adjacent 
sites on the surface is Oj = z/ ■ exp[—Eo{i)/kBT], where -E'o(^) is the activation energy for 
hopping of Xi atoms (or molecules). Here we ass ume that diffusion occurs only by th ermal 



hopping, in agreement with experimental results (IKatz et al. 19991 : iPerets et al. 20051 ). For 
small grains it is convenient to define the scanning rate, Ai = ai/S, which is approximately 
the inverse of the time it takes an Xj atom to scan the surface of the entire grain. 

The master equation provides the time derivatives of the probabilities P{Ni, N2, N3) 
that on a random grain there will be Ni adsorbed atoms/molecules of the reactive specie Xj. 
It takes the form 



3 

P(Xi, X2, Xs) = J2^^ [^(••' P(Xi, X2, X3)] 

i=l 

3 

+ [(^^ + + X,P(Xi, X2, X3)] 

i=l 
2 

+ J2 M{N^ + 2)(X, + 1)P(.., X, + 2, ..) - X,(X, - l)P(Xi, X2, X3)] 

i=l 

+ {A, + A2) [(Xi + 1)(X2 + l)P(Xi + 1, X2 + 1, X3 - 1) - XiX2P(Xi, X2, X3)] 

+ {A, + A,)[{Ni + 1)(X3 + l)P(Xi + 1, X2, X3 + 1) - XiX3P(Xi, X2, X3)]. (1) 

The terms in the first sum describe the incoming fiux, where Pj (atoms s"^) is the fiux 
per grain of the specie Xj. The second sum describes the effect of desorption. The third 
sum describes the effect of diffusion mediated reactions between two atoms of the same 
specie and the last two terms account for reactions between different species. The rate of 
each reaction is proportional to the number of pairs of atoms/molecules of the two species 
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involved, and to the sum of their scanning rates. The moments of P{Ni, N2, N3) are given 
by {N^N^N^) = 'En,,N2,N3 ^i^2NiP{Ni, N2, N:i), where a, 6, c are integers. In particular, 
(Ni) is the average population size of the specie Xj on a grain. The production rate per 
grain, R{Xk) (molecules s~^), of Xk molecules produced by the reaction Xi + Xj Xk is 
given by R{Xk) = (A + Aj){Ni Nj), or by R{Xk) = Ai{Ni{Ni - 1)) in case that i = j. 

In numerical simulations the master equation must be truncated in order to keep the 
number of equations finite. This can be done by setting upper cutoffs Nf^^^ i = 1, . . . , J 
on the population sizes, where J is the number of reactive species. However, the number of 
coupled equations, Ne = Yli=iiN^'''' + 1). 

grows exponentially with the number of reactive 
species. This severely limits the applicability of the master equation to interstellar chemistry 



( IStantcheva fc Herbst 20031 ). To reduce the number of equations one tries to use the lowest 
possible cutoffs under the given conditions. In any case, to enable all reaction processes to 
take place, the cutoffs must satisfy N^^^ > 2 for species that form homonuclear diatomic 
molecules (H2,02, etc.) and A^™^^ > 1 for other species. 

The average population sizes of the reactive species and the reaction rates are com- 
pletely determined by all the first moments and selected second moments of the distribution 
P{Ni, N2, N3). Therefore, a closed set of equations for the time derivatives of these first and 
second moments could provide complete information on the population sizes and reaction 
rates. For the simple network considered here one needs equations for the time derivatives 
of the first moments (iVi), {N2) and {N^) and of the second moments {Nf), (A^|), {N1N2) 
and {N1N3) [nodes and edges, respectively, in the graph shown in Fig. [T]. Such equations 
are obtained by taking the time derivat ive of each moment and using Eq. ([1]) to express 



the time derivatives of the probabilities (ILipshtat fc Biham 20031 ). Here we show two of the 
resulting moment equations: 



"^^^'^ = F, + {2A,~Wi){Ni)-2Ai{N^)-{A, + A2){NiN2)-{A, + A3){NiN3) 



dt 



dt 

- {A, + A:,){NiN^,) + {A, + A2){{N^N2)-{NiN2)-{N,N2Ns)). (2) 

In these equations, the time derivative of each moment is expressed as a linear combination 
of several other moments. However, the right hand sides of these equations include third 
order moments for which we have no equations. In order to close the set of moment equa- 
tions we must express th e third order moments in terms of first and second order moments 



(ILipshtat fc Biham 20031 ). This can be done by imposing the following constraint on the 



master equation: at any given time, at most two atoms or molecules can be adsorbed simul- 
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taneously on the surface. Furthermore, these two atoms or molecules must be from species 
that react with each other. The resulting cutoffs allow only eight non-vanishing proba- 
bilities, namely, P(0,0,0), P(0,0,1), P(0,1,0), P(1,0,0), P(2,0,0), P(0,2,0), P(1,1,0) 
and P(1,0,1). The third moments in Eq. can now be expressed in terms of these 
non-vanishing probabilities, giving rise to the following rules: (a) {N1N2N3) = 0; (b) 
(NfNj) = {NiNj); (c) {Nf) = 3{Nf) - 2{Ni). Using these rules, which are general and 
apply to any network of binary reactions, one can modify Eqs. ([2]), and obtain a closed set 
of the form 



d{Ni) 
dt 

dt 

dm 

dt 

d{NiNj) 
Jt 



F^ - W,{N,) - 2A{{N^) - (AT,)) - {A, + A,){N,N2) - k^Mi + 

P3 - W^3(iV3) - (^1 + A3)(iViiV3) + (^1 + A2)(iViiV2) 
+ (2P, + W^, + 4A,)(iV,) - {mi + 4A,)(iV5 

(Ai + A2)(iViiV2) - + A3)(iV,iV3) 

Pi(iV,) + P,.(iVi) - {Wr + W, + + A,)(iViiV,), (3) 



where i = 1, 2, j = 2, 3 and bij = 1 if z = j and otherwise. This set includes one equation 
that accounts for the population size of each reactive species and one equation that accounts 
for the rate of each reaction. Although these equations were derived using strict cutoffs, 
that are expected to apply only in the limit of very small grains and low flux, they provide 
accurate results for a very broad range of conditions. The point is that once the set of 
moment equations is derived, the probabilities do not appear anymore, so the constraint 
is not explicitly enforced. In fact, the equations maintain their accuracy even when the 
populations sizes of the reactive species are well beyond the constraints imposed above. 

In Fig. [11(a) we present the population sizes of H (+) and O (x) atoms on a grain, vs. 
grain diameter, obtained from the moment equations. In Fig. [UJ^b) we present the production 
rates of H2 (squares), O2 (triangles) and H2O (circles) vs. grain diameter, obtained from 
the moment equations. The results are in excellent agreement with the master equation 
(solid lines). In the limit of large grains they also coincide with the rate equations (dashed 
lines). Note that the moment equations apply even when there are as many as 10 hydrogen 
atoms on a grain. The parameters used in the simulations are s = 5 x 10^^ (sites cm~^). 
Pi = 5.0 X 10~^°S' (atoms s~^), P2 = O.lPi and P3 = 0. The activation energies for diffusion 
and desorption were taken as Po(l) = 44, Pi(l) = 52, Po(2) = 47, Pi (2) = 54, Po(3) = 47 
and Pi (3) = 54 meV. The grain temperature was T = 15K. The parameters used fo r 



hydrogen are the experimental results for low density amorphous ice (jPerets et al. 20051 ). 
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For the other species there are no concrete experimental results, and the values reflect the 
tendency of heavier species to bind more strongly. The fluxes and grain temperatures are 
suitable for dense molecular clouds. 



3. The Methanol Network 

Consider the case in which a flux of CO molecules is added to the network. This gives rise 
to the net work shown in Fig. [2l w hich includes the following sequence of hydrogen addition 



reactions (IStantcheva et al. 20021 ) : H + CO ^ HCO, H + HCO ^ H2CO, H + H2CO ^ 
H3CO and H + H3CO —>■ CH3OH. Two other reactions that involve oxygen atoms also take 
place: O + CO CO2 and O + HCO — > CO2 + H. This network was studied before using 
the multiplane method, which required about a thousand e quations compared to ab out a 



million equations in the master equation with similar cutoffs (ILipshtat fc Biham 20041 ). The 
moment equations include one equation for each node and one for each edge, namely, the 
network shown in Fig. [2]requires only 17 equations. We have performed extensive simulations 
of this network using the moment equations and found that they are in excellent agreement 
with the master equation. In Fig. [2](a) we present the moment-equation results for the 
population sizes of H, O and CO on a grain vs. grain diameter for the methanol network. 
In Fig. [2](b) we present the moment-equation results for the production rates per grain of 
some of the final products of the network. The results are in excellent agreement with the 
master equation (solid lines), and coincide with the rate equations (dashed lines) for large 
grains. The activation energies for diffusion and desorption of CO (Xj), HCO (Xg), H2CO 
(Xg) and H3CO (Xio) were taken as Eo{7) = 50, ^i(7) = 55, Eo{8) = 52, Ei(8) = 58, 
Eo{9) = 53, El (9) = 59, Eo{10) = 55 and Ei(lO) = 62 meV. The flux of CO molecules was 
taken as Fj = 0.2Fi. Note that experiments on CO desorp tion from ice surfaces indicate that 



-El (7) should be slightly higher than the value used here (ICoUings et al. 2004 ). However, it 
turns out that using a higher value would compromise the feasibility of the master equation 
simulations, which are required here in order to establish the validity of the moment equation 
method. 

When simulating highly complex networks one encounters the problem of obtaining the 
equations themselves. The ability to automate the construction of the equations becomes 
important. A crucial advantage of the moment equations is that they can be easily obtained 
using a diagrammatic approach. A detailed presentation of the diagrammatic method will 
be given in Barzel & Biham (2007). 
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4. Summciry and Discussion 

In summary, we have introduced a method, based on moment equations, for the simu- 
lation of chemical networks taking place on dust-grain surfaces in interstellar clouds. The 
method provides highly efficient simulations of complex reaction networks under the extreme 
conditions of low gas density and sub-micron grain sizes, in which the reaction rates are dom- 
inated by fluctuations and stochastic simulations are required. The number of equations is 
reduced to one equation for each reactive specie and one equation for each reaction, which is 
the lowest possible number for such networks. This method enables us to efficiently simulate 
networks of any required complexity without compromising the accuracy. It thus becomes 
possible to incorporate the complete network of surface reactions into gas-grain models of 
interstellar chemistry. To fully utilize the potential of this method, further laboratory exper- 
iments are needed, that will provide the activation energy barriers for diffusion, desorption 
and reaction processes not only for hydrogen but for all the molecules involved in these 
networks. 

We thank A. Lipshtat for helpful discussions. This work was supported by the Israel 
Science Foundation and the Adler Foundation for Space Research. 
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Fig. 1. — (a) The population sizes of H (+) and O (x) atoms on a grain vs. grain diameter, 
d (and the number of adsorption sites, S), obtained from the moment equations for the 
network shown above, where nodes represent reactive species, edges represent reactions and 
the products are indicated near the edges; (b) The production rates of H2 (squares), O2 
(triangles) and H2O (circles) on a grain vs. grain diameter, obtained from the moment 
equations. The results are in excellent agreement with the master equation (solid line). In 
the limit of large grains they also coincide with the rate equations (dashed line). 
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Fig. 2. — (a) The population sizes of H (+), O (*) and CO (x) for the methanol network 
shown above, vs. grain diameter, d (and the number of adsorption sites, S"), obtained from 
the moment equations, (b) The production rates of several final products of the methanol 
network per grain vs. grain diameter, obtained from the moment equations. The results 
are in excellent agreement with the master equation (solid line) . In the limit of large grains 
they also coincide with the rate equations (dashed line). 



